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ABSTRACT 


David S. Thoopaon, Haster of Science, 1980 

Hajor: Aerospace Engineering, Department of Aerospace Engineering 

Title of Thesis: Numerical Solution of the Navler-Stokes Equations 
for High Reynolds Number Incompressible Turbulent 
Flow 

Directed by: Dr. Joe F. Thompson 

Pages In Thesis: 53 Words In Abstract: 198 

Abstract 

The full Navler-Stokes equations for Incompressible turbulent 
flow must be solved to accurately represent all flow phenosiena which 
occur In a high Reynolds number Inexpressible flow. A two-layer 
algebraic eddy viscosity turbulence model Is used to represent the 
Reynolds stress terms In the time-averaged Incompressible Navler- 
Stokes equations in the primitive variable formulation. The development 
of the boundary-fitted coordinate systems has made the numerical 
solution of these equations feasible for arbitrarily shaped bodies. 

The non-dimensional time-averaged Navler-Stokes equations, 
including the turbulence model, are represented by finite difference 
approximations in the transformed (C,n) plane. The resulting coupled 
system of nonlinear algebraic equations is solved using a point 
successive over-relaxation (SOR) iteration. 

The test case considered in this study was an NACA 64A010 air- 
foil section at an angle of attack of two degrees and a Reynolds 
number of 2,000,000. Several boundary-fitted coordinate systems 
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were generated and uaed to evaluate various filters and various 
representations of the convective terms. Pressure distributions are 
presented which emphasize the difficulties associated with each 
technique. 

The preliminary results of a solution are presented which 
encourage the continuation of the solution to obtain a steady state 
solution. The major results of the evaluation of the techniques are 
also summarized. 
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I. INTRODUCTIOH 


Th« problem of accuretoly predicting the flotrfleld about an arbi- 
trary configuration In a high Reynolda numbar Incot^reaalbla flow haa 
provided aeemlngly Ineurmountable problema to reaearchera In thla area. 
The only way to represent fully all flow phenometta which occur at these 
conditions Is to solve the full Navler-Stokes equations for Incom- 
pressible turbulent flow. The primitive variable formulation must be 
used If multiple bodies or three-dimensional flow Is to be considered. 

Since the solution of the Navler-Stokes equations Is essentially a 
very coaq>lex boundary value problem, the validity of tha solution la 
depandant on the accuracy of the representation of the boundary values. 
If a conventional grid system is used for an arbitrary configuration. 
Interpolation will be required at the boundaries. Thla may lead to 
poor application of the boundary conditions. In high Reynolds number 
flow, there are large gradients in regions near solid boundaries. 

These gradients are generally dominant in determining the character of 
the solution. The resolution of these gradients requires that a large 
nund>er of closely spaced coordinate lines exist In the regions near 
solid boundaries. This would suggest using a fine mesh near these 
boundaries and a coarse mesh in the regions ^ere the gradients are 
small with some type of transitional mesh In between. 

A technique has been developed by Thompson, Thames, aiul Mastln 
[1] which will help to alleviate these problems. This technique 
numerically generates a discrete mesh system, called a boundary-fitted 
coordinate system, for arbitrary configurations. These mesh systems 
possess b. constant coordinate line coincident with each physical 


boundary so Intarpolatlon is alimlnatad at cha botmdarlas. By aodifying 
tha govarnlng aquations, coordinata linas can ba concantratad In any 
ragion of tha field. 

Tha two<>diBansional Maviar-Stokas aquations for incoaiprassibla 
turbulent flow ara reprasented by a finite diffaranca approxisMtion for 
tha tina-avaragad incomprassibla Naviar-Stokas aquations and a sli^tly 
modified version of tha algebraic turbulence model developed by Baldwin 
and Lomax [2]. Tha finite diffaranca approximations mist ba augmented 
by tha inclusion of terms relating tha discrata mash and tha physical 
grid. This effectively removes the physical coordinate system fr«a tha 
problm at the expense of complicating the original set of equations. 
However, application of the boundary conditions is simplified since 
the boundary conditions arc given on straight boundaries in the trans- 
forsMd plane. Since the finite difference approximations represent an 
elliptic system of nonlinear partial differential equations, an 
Iterativa technique must be used to obtain a solution. Tha technique 
used in this study was an accelerated Gauss-Seldcl iteration* or 
successive over-relaxation (SOR) . 

These techniques have been used with some success for incoaqiressible 
viscous flow by several researchers. Bearden [3] obtained results for 
laminar flow at a Reynolds number of 1,000,000 about a single element 
airfoil at zero angle of attack using the stream function-vortlcity 
formulation of the Navier-Stokes equations. Reddy [4] also obtained 
results for the same flow conditions uslne. the Integro-differential 
formulation. Primitive variable formu-.jtions have been developed by 
Hodge [5] and Shanks [6]. Hodge [5] obtained results for laminar flow 
about a single element airfoil at zero angle of attack for a Reynolds 
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nuabcr of 41,400. Shonko [6] considered low Reynolds nuiAer flow sbouc 
s substerged hydrofoil. 
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II. THE BOUNDARY-FITTED COORDINATE SYSTEM 


Much research has been devoted to the development of the techniques 
necessary for numerically generating boundary-fitted coordinate systems. 
Since the mathematical development and numerical Implementation of these 
techniques Is given In great detail by Thompson, Thames, and Mastln [1], 
Hodge [5], Thompson [7], and Thames [8], only an overview will be 
presented here. In addition, a method used to contract coordinate lines 
near a body in the field Is given In Appendix B. 

Consider a two-dimensional doubly-connected region as shown In 
Figure 1. The general transformation is one which associates each 
point (x,y) in the physical plane with a corresponding point (Ct<l) In 
the transformed plane. Let n “ n, on “he body contour and n ■ H2 
on the outer boundary The contour in the physical plane maps 

to the contour F* in the transformed plane. Similarly, the contour T 2 
maps to the contour F*. The contours F^ and F^ represent a "cut" to 
be made In the physical plane and constitute the reentrant segments, 

F^ and F^, in the transformed plane. Let 5 « on F^ and 5 ■ ^2 

on Fj. C Is allowed to vary monotonically from to ^2 both the 
inner and outer boundaries, F^ and F 2 respectively. The values of the 
physical coordinates on F^ and F^ are the same, but the function 
5 = C(x,y) is multivalued on F^ and F^ since ^ ? 2 * 

Now C and n have been completely specified on all the boundaries 
of a closed field. It remains to define the values in the interior of 
the field in terms of these boundary values. This implies that elliptic 
partial differential equations can be used to generate the field points 
since the solution of an elliptic partial differential equation is 
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r 




completely defined in the interior of a region by its values on the 
boundaries of that region. The elliptic system chosen must exhibit 
certain ^lav ^^nlllll principles which preclude the occurrence of extrema 
in the interior of the region. This will assure that a one-to-one 
correspondence exists between the physical and the transformed plane. 

The generating system of equations used in this study is a slightly 
modified version of the elliptic systems given by References [1], [5], 
[7], and [8]. The elliptic system used to generate the boundary- 
fitted coordinate system is given by 


n + 

XX yy 


. JL 


Q(C,n), 


(2.1b) 


subject to the following Dirichlet boundary conditions 


5 


Cj^(x.y) 

n 


'^l 

C 

B 


n 

- - 


"2 


[x.yJeFj^ 


[x.yjeF^ 


(2.2a) 


(2.2b) 


where P(5,n) and Q(5,n) are the attraction functions for the 5 and n 
lines respectively and a, y and J are given, along with other quantities 
relating the physical and the transformed planes in Appendix A. Since 
it is desired that all numerical computations be performed in the 
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transformed plane, the Independent and dependent variables must be 
interchanged. In the transformed plane, the generating system Is 
given by 

- -(ax^P(^;,n) + YX^^Q(5,n)) (2.3a) 

- 26y^^ + ■ -(ay^P(C,n) + Yy^Q(5,n)) (2.3b) 

with the transformed boundary conditions 


X 

at 

fl( 5 ,n^) 

. [ 5 ,nil£rJ 

(2.4a) 

y 


^2(5,113^) 



X 

m 

81(5,12) 

. [ 5 ,n 2 ]er* 

(2.4b) 

y 


§2 ^5 1 ^2 ) 




where the definition of 6 is given in Appendix A. The functions 
fj^(£,nj^), f 2 (C,nj^), gj^(C,n 2 ). and § 2 ( 5 , 02 ) are specified by the known 
shape of the contours and respectively, and the specified 5 
distribution thereon. 

Even though this system of quasi-linear partial differential 
equations, Equations (2.3a) and (2.3b), is more complicated than the 
original system, Equations (2.1a) and (2.1b), the boundary conditions 
are specified along straight boundaries in the transformed plane. 

Also, the coordinate line spacing in the transformed plane is uniform. 
At this point it should be noted that the actual values of 5 and n are 
irrelevant. The only quantities required by the finite difference 
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expressions are the values of AC and An which are taken to be unity by 
construction since cancellation occurs upon substitution Into the 
finite difference expressions. 

The generating system of equations. Equations (2.3a) and (2.3b), 

Is solved In the transformed plane. All derivatives are approximated 
by second-order central finite difference expressions. The resulting 
set of nonlinear simultaneous difference equations Is solved using a 
point SOR Iteration. 

Due to the Instability In the Navler-Stokes solution near the 
trailing edge reported by Steger and Bailey [9] and Thompson [10] 
for 0-type coordinate systems, the coordinate systems used In this 
study were generated using a different outer boundary configuration 
than shown in Figure 1. All coordinate systems used In this study 
possessed a "C-shaped" outer boundary as shown in Figure 2. The use 
of this configuration eliminates the problem of the coordinate lines 
having to "bend" around the sharp trailing edge. 

Consider the two-dimensional doubly-connected region shown In 
Figure 2. Once again the body is represented by the closed contour 
However, the "C-shaped" outer boundary is represented by three contours, 
and the downstream boundaries and F^. The cut in the physical 

plane is made along the contours F^ and F^. For this configuration, 

n = n, on the contours F^, F, and F^ and n = Ho on the contour r». £ 

1 D 1 o z 3 

varies monotonically from £ onr 2 to£ ~ ^2 again 

the contour in the physical plane maps to the contour Fj^ in the 
transformed plane, 1'2 maps to ®tc. Since the value of n is constant 
along the contours F^, and F^, these contours must represent a line 


7 


of constant n In the transformed plane. Also, the cut made in the 
physical plane along the contours and F^ is represented by the re~ 
entrant segments F* and F* in the transformed plane. 

Several different forms of the attraction functions, P(€»n) and 
Q(C>n) from Equations (2.3a) and (2.3b), are given by References [1], 

[3], [4], and [7]. In this study only n~Hne contraction was used, so 
P(^,n) Is taken to be zero. Initially the form of the attraction 
function, Q(C,n)» was taken to be the form given In the appendix of 
Reference [2]. The form of Q(C,n) which gave the best results, however. 
Is presented In the appendix of Reference [4] and Reference [7]. This 
technique Is developed in detail in Appendix B. 
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111. THE NAVIER-STOKES EQUATIONS IN PRIMITIVE VARIABLES 

The tlae-averaged Navier-Stokes equationa for Kwo-dinenalonal In- 
coflq>resslble flow coupled with an algebraic eddy viscosity turbulence 
model are presented as an alternative to the full Navier-Stokes 
equations for two-din^nsional incompressible turbulent flow. The 
primitive variable formulation is employed. The resultant equations 
are non-dimensionallzed and transformed to the transformed (C«n) plane. 
The boundary conditions and their application are also presented. 


A. The Basic Equations 

As stated previously, the full Navier-Stokes equations for in- 
compressible turbulent flow must be solved to accurately predict a high 
Reynolds number incompressible flowfleld. Since an extension to 
multiple bodies and eventually three-dimensional flow is desired, only 
the primitive variable formulation will be considered. The full 
Navier-Stokes equations for two-dimensional incompressible flow are 
given in the primitive variable formulation by 


+ 

uu + 

vu ) ■ 

- p 

+ 

P (u 


u 


X 

y 

• X 

XX 


yy 

+ 

uv + 

w ) ■ 

- D 

+ 

p(v 


V 


X 

y 

y 

XX 


yy 



D u 

+ V ■ 

0 




(3.1a) 
(3.1b) 

X y <3.1c) 

where u and v are the velocity components parallel to the x and y 
directions respectively, p is the pressure, p is the fluid density, v 
is the molecular viscosity coefficient, and the subscripts x, y, and t 
represent partial differentiation in the usual manner. Equations (3.1) 
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theoretically Include the turbulent motion of the fluid If the time 
Btep size and the spacing of the discrete mesh points are taken to be 
arbitrarily small. This approach Is Impractical due to the excessive 
conq>utatlonal requirements. Some approximate method must be used to 
model the effects of turbulence. 

The time-averaged Navler-Stokes equations for Incompressible flow, 
given In the primitive variables by 


p(u^ + u5 + ^„) - -p„ 

+ 


-p(u'v') -p(u^) 

(3.2a) 

t X y X 

XX 

yy 

y X 


+ 


-p(u'v') -p(v'^) 

(3.2b) 

p(v„ + uv„ + vv. ,) ■ -p„ 

+ V ) 

t X y y 

XX 

yy 

X y 

D - 

U + V 
X y 

- 0 


(3.2c) 


where the over bars Indicate time-averaged quantities and the primes 
Indicate fluctuating quantities, were considered In this study. This 
form was chosen because of the availability of techniques which model 
the Reynolds stresses, -p(u'v'), -p(u'^), and -p(v'^) . 

The major problem associated with the primitive variable formu- 
lation of the Navler-Stokes equations. Equations (3.1), Is the lack of 
a time derivative for pressure. There Is no direct way of advancing 
pressure to the next time level. In fact, the role of pressure In 
incompressible flow Is to somehow adjust Itself so that continuity 
will be satisfied. A Poisson equation in pressure can be obtained by 
taking the divergence of the momentum equations. Equations (3.1a) and 
(3.1b). Several forms of this equation are given by Hodge [5] and 
Shanks [6]. The form used In this study Is given by 

"(P + p ) * (u )^ + 2v u + (v )^ + D (3.3) 

*^xx *^yy X X y y t 
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where Is the tine derivative of the continuity aquation. Analyti- 
cally* D la always zero, so will always be zero. Numerically, this 
is generally not true. The term is retained In Equation (3.3) as a 
corrective term to adjust the pressure in an attempt to drive 
continuity to zero. 

B. The Turbulence Model 

The technique used in this study to model the Reynolds stresses 
is a slightly modified version of the two-layer rlgebralc eddy viscosity 
turbulence model developed by Baldwin and Lomax [2]. In this model, 
the effects of turbulence are represented by an eddy viscosity co- 
efficient That is, the Reynolds stress terms of Equations (3.2a) 
and (3.2b) are dropped and the molecular viscosity coefficient p is 
replaced by p + p^. Equations (3.2) then become 

p(u^ + uu + vu ) ■ -p + p(l + c)(u + u ) (3.4a) 

t X y *^x XX yy 

o(v. + uv + w ) ■ -p + p(l + e)(v + V ,) (3.4b) 

t X y y XX yy 

D ■ u + V “0 (3.4c) 

X y 

where c is the ratio of to p. In addition, the distribution of 
vortlclty is used to determine the length scale so that finding the 
edge of the boundary layer is not necessary. 

Spatial derivatives of the eddy viscosity have been neglected 
in both the momentum equations and in the Poisson equation for the 
pressure. This approximation was applied in order to avoid consider- 
able complication of these equations, but justification was not 
established. 
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The method ueed to model treneitlon Is given by Cebeci end 
Bradshaw [11]. The calculation of the eddy viscosity coefficient la 
modified by an Intermittency factor that accounts for the transitional 
region that exists between laminar and turbulent regions of a flow. 

This avoids the assiunptlon that the laminar flow becomes turbulent at 
the transition point which can in general lead to substantial error. 

For simplicity of calculation, the body is assumed to be a flat plate. 

A transition point is calculated on the upper and lower surfaces by 
assiunlng that transition occurs at the first minimum pressure point 
on each surface. 

C. The Non-Dimensional Equations in the Transformed Plane 

Equations (3.4) and Equation (3.3) can be simplified considerably 
by using the following non-dimensional variables: 

X* ; x/L, y* = y/L, 

u* = u/U , V* = v/U_,, 

OD ww 

p* = (i - pJ/pU^„, t* E tu„/L, 

Re i pU^L/w 

where U is the freestream velocity, L is the characteristic length, 

00 

p is the freestream pressure, and Re is the Reynolds number. After 
substitution of the non-dimensional variables, with the asterisks 
dropped for convenience henceforth. Equations (3.4) and Equation (3.3) 
become 

u + uu + vu “ -p + (u + u ) (1 + c)/Re (3.5a) 

t X y *^x XX yy 

V + uv + w ■ -p + (v + V ) (1 + e)/Re (3.5b) 

t X y y XX yy 
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D ■ u + V ■ 0 
X y 


(3.5c) 


“(p + P ) ■ (u )^ + 2v u + (v + D^. (3.5d) 

'*^xx *^yy X X y ' y t 

V 

This set of equations represents the form used for the bulk of this 
study. 

Equations (3.5) must now be transformed Into the transformed plane 
using the relations and definitions given In Appendix A. The resulting 
transformed equations are given by 

"t ■ 

■ * “““« ■ “"tn * 

+ (o/J^)u + (t/J^)u- 1 (1 + e)/Re (3.6a) 

n t 


't • hV'-’ * ’'*’‘5% ■ 




+ (o/J^)v + (t/J^)v ] (1 + e)/Re 
n ^ 


D - (y^uj - y^u^ + X|.»^ - - 0 


(3.6b) 

(3.6c) 


-[(ap^^ - 20p^^ + + (o/j2)p^ + (t/j2)P^J 


■ Vn><Vn ■ V5> 


+ (x V - X v-)2]/J^ + D 
4 n n ^ t 


(3.6d) 


Equations (3.6) are now given on a rectangular field with a square grid 
In the transformed plane. The numerical procedures used to obtain 
solutions to Equations (3.6) are given In Chapter 4. 
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D> The Boundary Conditions 


The boundary conditions on the surface of the body are given by 

u - 0 (3.7a) 

V - 0 (3.7b) 

which are the no~sllp boundary conditions for a vlscou*t fluid on the 
surface of a stationary body with no transpiration. However, the 
pressure on the wall p^ Is unknown and must be calculated. Hodge 
[5] presents two methods for obtaining the pressure on the body surface. 

A Chorln-type pressure Iteration utilizing the continuity equation 
to obtain the wall pressure Is given by 


(s+l) , (a) _ 


QD 


<3. 8a) 


"w "w 

where s Is the Iteration counter and il Is an appropriate acceleration 
parameter given by Hodge [5] to be 


n - 0ij2/t2At(a+Y)] (3.8b) 

where u is an acceleration parameter and At Is the time step size. 

Another method of obtaining the wall pressure Is to evaluate the 
normal derivative of the pressure at the wall from the momentum equation. 
The pressure normal to an n-wall is given by 

• Vp • (yp_ - 6p.)/J»'7 (3,9) 

n n ^ 

where n^^^ • Vp is found from Equation (A. 27) in Appendix A. The 
normal component of the momentum equation for a body with no trans- 
piration and the nc-slip boundary condition is given by 

- [BP^ + /^RX - x^RY]/y (3.10a) 
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wh«r« 


RX s -[(-23u. + Yu„„)/J + (o/J)u„]/R« (3.10b) 

RY 5 -((-2ev.„ + Yv^„)/J + (o/J)v„l/R* (3.10c) 

tn nn n 

Th« wall pressure can be calculated from Equations (3.10) using one- 
sided difference equations. 

The pressure at a sharp trailing edge la calculated by applying 
Equations (3.8) or Equations (3.10) to the upper and lower surfaces of 
the body at the trailing edge. Since these two values are not generally 
equal, an average la taken tr' obtain the pressure at the trailing 
edge. 

The frees tream boundary conditions are applied at all points on 
the outer boundary. Ttie freestream conditions are given In terms of 
the non-ditncnslonal vrjrlables by 


u ■ cos 


(3.11a) 


V ■ sin 


(3.11b) 


p - 0 


(3.11c) 


where ^ is the angle of attack. The freestream boundary conditions can 
be appX'‘ed on the downstream boundary because of the large distance from 


the trailing edge to the downstream boundary used in this study. 


IV. IMPLEKENTAT1(»I OF THE NUMERICAL SOLUTXOf 


Th« transformed equations must now be approximated by finite 
difference expresaions. An outline of the technlquea uaed to develop 
a finite difference approximation is preaentad. The general forsw of 
the finite difference expressions used in the approximation arc given 
in Appendix C. The method of solution of the finite difference 
approximation is discussed. The finite difference approxiauttion 

for derivatives across the cut is discussed. Various numerical 
techniques used to improve the stability of the numerical solution are 

diacussed. 

A. The Finite Difference Approximation 

The transformed Navier-Stokes and Poisson pressure equations, 
Equations (3.6), must now be represented by a finite difference 
approximation on the discrete mesh system. Since this finite difference 
approximation represents a system of nonlinear partial differential 
equations which are elliptic in space and parabolic in time, the 
finite difference approximation choben must accurately reflect these 
characteristics. Also, consideration must be given to the stability 
requirements of the method used. 

For these reasons, an Implicit algorithm was developed which 
utilized backward-time and cen ro -space differencing techniques. 

Both first-and second-order differences were considered for the time 
derivatives. Only second-order differences were UoCd for the spatial 
derivative!. Since central differences cannot be used to represent 
spatial derivatives in the wall pressure calculc’ ion, Lquatlon (3.8a) 
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or Equation (3.9). aacond-ordar forward dlffarancas vara uaad to 
rapraaant cha n-darivacivaa. Tha ganaral foras of cha finlca 
diffaranca axprasalons uaad in thia algorithm ara givan in Appandix C. 

At thia point, it ahould ba noted that darivativaa taken acroaa 
tha "cut" in tha phyaical plane muat receive apecial treatment. With 
rafaranca to Figure 2, a derivative taken acroaa tha contoura and 
in tha phyaical plana ia taken acroaa tha reentrart aagmanta r* and r* 
in tha tranaformad plana. Since x and y ara equal along thaaa linaa, 
evaluation ia naceaaary along only one of them. Aa ahown in Figure 3, 
pointa in the tranafomed (C,n) plane which are located at (Il-N,l) 
and (I24N,1) for N ^ II have the aame coordinataa in the phyaical 
(x,y) plane. With thia in mind, the finite difference uxpreaaiona 
given in Appendix C can be used to obtain 


^^n^Il“N,l " ^^n^I2+N,l * ^^Il-N,2 " ^I2+N,2^^^ (^-l*) 

^^m^ll-N,l ■ ^^nn^i2+N,l * ^li-N,2 " ^^il-N,l ^ll+N,2 


^^Cn^Il-N.l “ ^^Cn^I2+N,l 


^^Il-N-1,2 ■ ^Il-N+1,2 ^I2+N+l,2 


■^I2+N-l,2^''^ 


(4.1c) 


where f is an arbitrary function of C and n. Th'.e insures that the 
derivatives taken across the reentrant segments rre continuous. 

The nonlinear system of algebraic equations which is formed by 
the application of the appropriate finite difference expressions from 
Appendix C to Equations (3.6) is solved using a point SOR iteration. 
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The calculation of a linearly optimum acceleration parameter for an 
Iteration of this type followi! the procedure given by Thon^son (7]. 


B. Numerical Stabilizing Techniques 

Several purely numerical techniques were Implemented In an attempt 
to obtain a stable solution. These techniques are needed because of 
the nonphysical oscillations of the dependent variables which can 
develop In regions of steep gradients. These oscillations, or 
numerical Instabilities, can cause an otherwise stable solution to 
become unstable and diverge. 

There are two possible approaches which can be taken to enhance 
the stability of a numerical solution. One possible approach Is to 
treat the symptoms of the nrmerlcal Instability by filtering or 
equivalently, adding an artificial viscosity. The other approach Is to 
attempt to alleviate the cause of the Instability. 

A basic filter considered In this study Is the switched form of 
the Shuman filter given by Harten and Zwas [12]. The two-dimensional 
Shuman filter for a general function f Is given by 


^ “ f "b “ f f 

•i.J ^i,j ^ 8 ^''i+l.j 




where 6 » 1 for the simple Shuman filter. Since this technique has the 
effect of adding an artificial viscosity, the effective Reynolds 
number will be lowered considerably if the Shuman filter is applied 
at all points in the field. However, this will eliminate all of the 
oscillations. The characteristics of this filter can be Improved 
if it is locally applied only when a certain waveform is encountered. 
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A 


V 


If an N-waveform is to be filtered, 9 would be set to unity if an 
N-waveform Is encountered but would be set to zero otherwise. Thus, 
only the two sharp points of the N-waveform would be filtered. 

Similar statements can be made In relation to a filter for W-wave- 
forms or a W-fllter. Only the central point of a W-waveform la 
filtered. It should be noted that the N-fllter, and especially the 
U-fllter, add no diffusion at all to most of the field. The purpose 
of the N-fllter and the W-fllter Is not to eliminate the oscillations, 
but to control them so the solution does not diverge. 

The filter is applied after each time step. It could be applied 
after each iteration but this often causes convergence difficulties. 

In addition, the technique of filtering can be applied repetitively 
after each time step to further reduce the amplitude of the oscillation. 

Another technique which introduces an artificial viscosity is a 
fourth-order smoother related to that used by Baldwin and MacCormack 
[13]. Several forms of this smoother were implemented. The form which 
was finally considered simply replaces (1+c) in Equations (3.6) with 
(1+e+V^pJAtRe) , where V^p is the Laplacian of the pressure and At is 
the time step size. This technique introduces a significant artificial 
viscosity only in regions where the Laplacian of the pressure is large. 

Most of the attempts to alleviate the cause of the instability 
center around the representation of the convective derivatives. The 
form given by Equations (3.1) is the non-conservative form. The term 
uu^ + vu^ in Equation (3.1a) is actually the expanded form of 
(u‘)^ + (uv)y after cancellation of the continuity equation, 

Equation (3.1c), with D = 0. Analytically these two terms are 
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equivalent, numerically they may not be. Four different techniques 
were tested and are outlined here. Only terms In the x-im>mentum 
equations are shown since extension to the y-momentum equation Is 
similar. 

Two possible representations of the convective derivatives are 
noted In MacCormack [14]. The convective derivative can be 

represented as the average of a product (AOP) given by 


(u2)^ - i • (^-3) 

Equation (4.3) can itself cause Instabilities as explained In Reference 
[14]. A more stable technique Is the representation of the convective 
derivative as the product of an average (POA). This form Is 

given by 

^^'^i+1 

These techniques are described In greater detail by Reference [14]. 

An attempt to improve the stability of the techniques given by 
Equation (4.3) and Equation (4.4) was made by including a mass 
residual correction. This technique replaces uu^ + vu^ in Equation 
(3.1a) with 


(u^) + (uv) - u(u + V ) 


(4.5) 


X y X y 

which is analytically equivalent to the non- conservative form. The 
terms (uv)^ can be replaced by terms of the form of Equation 

(4.3) or Equation (4.4). 

Beam and Warming [15] used a technique to linearize the convective 
terms to prevent the Instability. The form used in this study is 
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given by 


- (u2)“ + 2u"(u"'^^ - u") (A. 6) 

where n Is the time step nusd>er. Equation (4.6) is a second-order 
Taylor series expansion for about (u^)**. Since u** Is known, 

this representation for (u^)“ is linear. 


21 


V. RESULTS 


The conputer code used in this study was written by Joe F. Thompson 
as part of the current research of the Department of Aerospace 
Engineering at Mississippi State University for the NASA Langley 
Research Center. The test case considered in this study was the flow- 
field generated by an NACA 64A010 airfoil section at an angle of 
attack of two degrees and a Reynolds number of 2,000,000. 

All the coordinate systems considered in this study were generated 
using the method of Chapter 2. The C-attractlon function P(C>n) was 
obtained by the simultaneous solution of Equations (2.3) after the 
initial guess had been obtained. The initial guess was fomred by 
placing the n-llne distribution produced by the contraction near the 
body on each 4~line. Now the ^-attraction function is given by the 
product of P(C,n) as obtained above and a decaying exponential based 
on "0.2d where d is the distance from the body. The n~attractlon 
function Q(C>n) was varied to obtain different n-line distributions 
near the airfoil. In all cases, the coordinate system was "C-shaped" 
with 113 C-lines and 51 n-lines in the field. There were 72 unique 
points on the airfoil. The downstream and outer boundaries were 
located 10 chord lengths from the airfoil. In addition, a Neumann 
boundary condition was Imposed on the n~lines at the downstream 
boundary. A typical coordinate system generated during this study is 
shown in Figure 4. The region of the field near the airfoil is shown 
In Figure 5. The convergence criteria used for all coordinate systems 
was 10“^. 
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Similarly, the following techniques were common to most of the 
Navler-Stokes solutions att^pted during this study. 

(1) A linear start in 100 steps was used with time step size 
of 0.01. 

(2) The first time step was run using first-order time differencing 
with no turbulence. 

(3) The basic turbulence model without the transition point 
calculation outlined in Chapter 3 was initiated along with 
second-order time differencing after the first time step. 

(4) The non-conservative forms of the momentum equations, 

Equations (3.6a) and (3.6b), were considered without the 
stabilizing techniques described in Chapter 4. 

(5) A zero first time derivative projection for the initial guess 
at the next time level was used. 

(6) The pressure on the boundary was obtained from Equations (3.8) 
with 01 >• 0.5. 

(7) The field pressure acceleration parameter was 1.0. The 
velocities were multiplied by 2.0 in the calculation of an 
optimum acceleration parameter. 

(8) The velocity and pressure convergence criteria were 10”^ and 
10““* respectively. 

(9) Partial convergence was accepted after 100 iterations. 

(10) All computations were performed on a UNIVAC 1100 series 

computer. 

Exceptions to these procedures are noted in the course of the discussion 
of results which follows. 
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The first type of coordinate system considered in this study was 
generated by forming a composite functloii for Q(C»n) es described by 
Thompson [7]. A Blaslus boundary layer profile was joined to a 
quartlc polynomial at the second line inside the boundary layer. The 
first coordinate system generated using this technique (CSl) had 10 
n-llnes contracted Into the boundary. The n'Hne distributions In 
the boundary layer for CSl at the point of maximum thickness of the 
airfoil and the leading edge of the airfoil are shown in Figures 6. a 
and 6.b respectively. 

CSl was used to generate a solution which diverged at t ■■ l.OS. 

The divergence occurred near the leading edge at approximately the 
first n-llne off the body. As shown In Figure 7, a pressure peak 
also occurred near this point. It was thought that the divergence 
was caused by an Insufficient number of n-llnes being contracted 
into the boundary layer. 

To remedy this problem, a second coordinate system (CS2) was 
generated using the same technique which had 20 n-llnes contracted 
into the boundary layer. CS2 was used extensively to test the various 
numerical techniques implemented during this study. The n-line 
distributions for CS2 are shown in Figures 8. a and 8.b. As shown 
in Figure 8.b, the n-lines are contracted much closer to the airfoil 
at the leading edge than previously shown in Figure 6.b . 

CS2 was used to generate a solution which diverged at t ■ 1.08. 

The divergence occurred in the region of the leading edge at the eighth 
n-line off the body. Once again, a pressure peak occurred in the 
region near the divergence as shown in Figure 9. This solution wac 
restarted from t ■ 1.00 using the W-f liter as described in Chapter 4. 
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Once again, the iteration diverged at t “ 1.08 and the pressure behavior 
was the same. Repeating this procedure using the N>filter described 
in Chapter 4 also had little effect on the divergence of the iteration 
or the pressure behavior. Also, the repetitive application of the 
W-fllter after each time step had little effect on the behavior of the 
solution. The solution was again restarted from t ■ 1.00 but with a 
new initial guess for each time step. The initial guess at the new 
time step was given by the previous time step solution. This also had 
little effect on the divergence or the pressure behavior. Applying the 
W-f liter from t ■ 0.00 also had no significant effect. Applying the 
simple Shuman filter. Equation (4.2), with 6 * 1, to the restart of this 
solution produced a more rapid divergence. 

An attempt was then made using CS2 to generate a solution using 
the momentum equation. Equations (3.10), to obtain the pressure on the 
boundary. The iteration diverged at t 1.08 and exhibited the same 
behavior as before. The pressures on the airfoil obtained using 
the momentum equation were essentially the same as those obtained using 
the continuity equation. Restarting the solution from t ■ 1.00 and 
using the N-filter had little effect on the divergence of the pressure 
behavior. Once again, taking the initial guess at the new time step 
to be the previous time step solution had no significant effect. 

Since each of these attempts was diverging soon after the start 
up procedure was completed, several different start-up procedures were 
tested using CS2 as the coordinate system. A 100 step cosine start was 
implemented. The iteration diverged at t ■ 0.86. The cosine start 
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VM no longer considered since It diverged before the lineer start. 
Another start-up procedure tested was the use of first-order beckvard 
tine differences for the tine derivatives. This also proved unsuc- 
cessful, as the iteration diverged at t ■ 0.94. This was sooner than 
the divergence of the second-order backward time differencing previously 
used. In addition, a 100 step linear start was used with the velocity 
multiple being 1.0 Instead of 2.0 in the optimum acceleration parameter 
calculation. This had a destabilizing effect on the iteration. 
Instabilities appeared In the solution as early as t “ 0.80. Divergence 
occurred at t - 0.08. The pressure behavior was essentially the same 
as stated previously. 

Several attempts were made to stabilize the solution by repre- 
senting the convective terms as described in Chapter 4 using CS2 as the 
coordinate system. The product of the average representation of the 
convective terms (POA) , given by Equation (4.4), was Implemented in 
an attempt to stabilize the solution. This Iteration diverged at 
t ■ 0.74. Convergence was lost as early as t <■ 0.50. Once again, 
the pressure behavior was similar to that shown In Figure 9. The mass 
residual correction, given by Equation (4.5), was used In conjunction 
with the POA representation of the convective terms In an attempt to 
stabilize the iteration. The combination of these techniques proved 
to be unsuccessful since this iteration diverged sooner than the 
iteration without the mass residual correction. The average of the 
product representation of the convective terms (AOP) , given by Equation 
(4.3), was implemented in conjunction with the mass residual correction. 
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given by Equation (4.5), In an attempt to atablliae the eolutlon. 

Thle wae alao unauccesaf ul , ae the Iteration diverged at t ■ 0.54. 

The AOP repreaentatlon waa not coneldered without the maaa reeldual 
correction becauae of the nonlinear instability reported by HacCormack 
[14]. 

In addition, a second-order linearisation of the form given by 
Equation (4.6) was implemented to enhance the stability of the solution. 
This was unsuccessful in that the divergence of the iteration occurred 
at t " 0.53. The divergence of the iteration once again occurred in 
the saiM region. 

Since the Iteration was diverging In the same region of the field 
for all cases using the coordinate system CS2, a coordinate system was 
generated (CS3) with 20 n-llnes contracted into the boundary layer 
using the attraction function given by Roberts [16]. This attre tion 
function is a continuous exponential function. ’H-.is eliminates any 
problems associated with discontinuities in the derivatives of the 
n-line distribution which result from a composite function being 
chosen for the attraction function. Figures 10. a and 10. b show the 
distribution of n-llnes in the boundary layer at the point of maxlimim 
thickness and the leading edge of the airfoil respectively. 

CS3 was used to generate a solution which did not diverge and 
seemed to be converging at time t * 2.00. However, this solution was 
worthless as shown In Figure 11. This Invalid pressure distribution 
was obtained from an equally Invalid velocity field. Thle behavior 
was caused by the looseness of the coordinate system convergence 
criteria. The distance at the leading edge to the first n-llne off the 
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c waa 2 x 10~^, which ia tha saaa order of magnitude as the con- 
vergence criteria. This problem could be alleviated by generating a 
coordinata system with a decrs(i*ed convergence criteria or by 
decreasing the nund>er of n-linea contracued into the boundary layer. 

The latter approach was implemented and a coordinate system (CS4) 
was generated with 10 n-lines contracted into the boundary layer. The 
distance at the leading edge to the first n-line off the body was 
2.6 X lO'**. The distributions of the n-lines in the boundary layer at 
the point of maximum thickness and leading edge of the airfoil are 
shown in Figures 12. a and 12. b respectively. 

CS4 was used to generate a solution which diverged at t ■ 1.18. 

This solution was more well-behaved than the solution obtained using 
CS3. The pressure distribution obtained at t ■ 1.00 is shown in 
Figure 13. Once again. a pressure peak occurred in the region where the 
divergence originated. Application of the W-fllter to the restart of 
this solution had no significant effect. 

A new coordinate system was generated (CS5) which utilized the 
techniques outlined in Appendix B. The distributions of the n-lines 
in the boundary layer at the point of maximum thickness and the 
leading edge of the airfoil are shown in Figures 14. a and 14. b 
respectively. As an attempt to reduce the roundoff error, o and t 
are calculated directly from P(C,n) and Q(C,n). This Iteration 
diverged at t “ 1.08 and exhibited the same behavior as observed for 
CS3. The MacCormack fourth-order smoother as described in Chapter 4 
was Implemented in an attempt to stabilize the solution. This Iteration 
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had mot divargad at t ■ 2.00, but eonvargmca waa loat aa aarly aa 
t *• 1.90. Tha aaxiana itaratlra arror occurrad Mar tha uppar aurfaca 
trailing adga at tha alxtaanth n**llna off tha body. Tha praaaura 
dlatrlbutlon obtalnad at t ■ 1.90 la ahoun In Flgura IS. Tha poor 
repraaantatlon of tha flovflald was cauaad by tha turbulanca nodal 
inaccurataly rapraaantlng tha addy vlacoalty. Tha addy vlacoalty 
coafflclant obtalnad from tha turbulMca nodal waa ao larga that ttM 
flow aaparatad at tha laadlng adga on both tha uppar and lowar aurfacaa. 
Ilila aaparatlon waa cauaad by tha aaau^tlon that tha antlra flowflald 
waa turbulant. Howavar, tha flow naax tha lairing adga waa actually 
lanlnar. Thla lad to tha Inclusion of tha transition point calculation 
aa diacussad in Chapter 3. 

At thla point in tins, adaquata conputatlonal facilltias bacana 
unavailable locally. This nacasaltated a c<xivaraion of all prograsM 
for conpatlblllcy with tha CDC computational facllltlaa available at 
NASA Langley Research Center. Because of this delay, only prallnlnary 
data is available from the new techniques Inplenantad during tha 
conversion. 

CSS was regenerated with a convergence criteria of lO'^. Using 
this new coordinate syatem, a solution was generated by including 
the transition point calculation outlined in Chapter 3, in addition to 
Che fourth-order MacCormack smoother and the calculation of o and t 
from and QU.n). Preliminary results indicate chat the transition 

point on the lower surface is calculated to be near the trailing edge 
of Che airfoil. Previously, all of the lower surface was turbulent. 

Now all of the lower surface is laminar. Tnis trend occurred from 
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C ■ 0.00 to t • 0.70, the last tine lavol lAora data la avallabla. 
Also, full convarganca has baan obtained up to this tlaM lavol. 


VI. COHCLUSZOI^ 


The d«v«IoiMMtC of the boundary-fitted coordinate aysteBa have 
■ada the nuaMrieal aolution of the Naviar-Stokaa aquationa for high 
Seynolda nuadiar incoopraaaibla turbulent flow faaaibla. Several 
nvoMrical tachniquaa were iaiplenented and evaluated during thia 
atudy. The following liat atanariaea the noat iaq>ortant reaulta of 
thia evaluation. 

(1) The non-conaervative formilation of the Naviar-Stokaa aquationa 
givea better reaulta than foranilationa involving AOP or POA 
repreaentationa of the convective tenia. In addition, the 
■aaa reaidual correction and the aecond-order lineariaation 
were alao ineffective. 

(2) The MacComack aaoother was the only filter which enhanced 
the stability of the solution significantly. The Shuoan 
filter, N-filter, and W-fllter did not significantly affect 
the stability of the solution. The ineffectiveness of these 
filters can be attributed to the larger time step size of 
lapllcit methods. The filter was not applied as often 

over the same time interval as it would have been for the 
explicit methods for which it was developed. 

(3) The pressure on the boundary can be obtained from the continuity 
equation via a Chorin-type Itt^ratlon or from the component 

of the momentum equation normal to the boundary. However, 
using the continuity equation seems to Improve the stability 
of the solution slightly. 
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(4) Prellffltnary data indicates that the addition of an inter- 


mittency factor at the transition point significantly iiq>roves 
the effects of the algebraic turbulence model used in this 
study. 

(5) Better results are obtained for coordinate systems which 
have a continuous distribution of coordinate lines with 
small convert^ence criteria (10~^). 

The prellmlnaxry . esults of the final solution described in 
Chapter 5 are the most encouraging to date. This solution should be 
continued in an attempt to obtain a valid steady-state solution. 
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RECTANGULAR TRANSFORMED PLANE 

Figure 1. Coordinate Transformation - 0-Type Outer Boundary 
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Figure 2. Coordinate Transformation - C-Type Outer Boundary 
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Figure 4. Typical Coordinate System-NACA 64A010 
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Figure 6a. Eta-Line Distri- 
bution at Point of Maximum 
Airfoil Thickness - CS2 
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Figure 6b. Eta-Line Distri- 
bution at Leading Edge - 
CSl 
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Figure 8a. Eta-Line Distri- 
bution at Point of Maximum 
Airfoil Thickness - CS2 


Figure 8b. Eta-Line Distri- 
bution at Leading Edge - 
CS2 
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Figure 10a. Eta-Line Distri- 
bution at Point of Maximum 
Airfoil Thickness - L'S3 


Figure 10b. Eta-Line Distri- 
bution at Leading Edge - 
CS3 
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Figure 12a. Eta-Line Dlstri' 
but ion at Point of Maximum 
Airfoil Thickness - CS4 


Figure 12b. Eta-Line 
buticn at Leading E 
CS4 
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Figure 14a. Eta-Line Distri- Figure 14b. Eta-Line Distri- 
bution at Point of Maximum bution at Leading Edge - CSS 
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APPENDIX A 


Various Relations in the Transformed Plane 

This Appendix contains many of the pertinent definitions and 
relations in the transformed (Cin) plane. Thompson, Thames, and 
Mastin [1] present a comprehensive set of relations and the notation 
used in Reference [1] is retained here. Similarly, the following 
function definitions are applicable: 

f(x,y,t) - A twice continuously differentiable scalar function of 
x,y and t 

F(x,y) " i Fj^(x,y) + j F 2 (x,y) - A continuously differentiable 

vector- valued function^ 1 and j 
are the conventional Cartesian 
coordinate unit vectors. 

It should be noted that all derivative transformations given here are 
in the geometrically non-conservative form. 

Definitions of the Transformation 



a c x^ + y^ 
n n 

(A.l) 


P = x^x + y^y 

(A. 2) 


Y ^ 

(A. 3) 

Dx 

= * '■\n 

fA.4) 

Dy 

= + yy^^ 

(A. 5) 


0 = (y^Dx - x^Dy)/J 

(A.6) 


T H (x Dy - y Dx)/J 
n n 

(A. 7) 


J - x^y - X y. 

4 n n-^4 

(A. 8) 
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Where J is the detemlnant of the Jacobian matrix 


Derivative Trane format ions 


XX 


'x ' ■ ‘Vs - 

*, = ■ ‘»s‘n ■ 

*<''n''S6 - * i^n>‘%'s ' ’‘s‘n>'-'' 

+<yJ«j5 - 2y5V^„ + y|*,„) 




(A. 11) 


yy 


(3 


Xjt 5 n % nn 

+(x^y__- 2x-x y. + xfy ) (x f. 

C t'nn n K 

■ ^“‘s-n-Sn * *Kn> 

<>s% - Vs>‘-'' 


x^t„)/j2 


(A. 12) 


‘xy • “Vn * *n^> Ur, - ViUn - Vn‘s5>"' 
‘%%’‘S6 ■ ‘’‘s'’n Vs^ Un * V(\r,’(UU 
* ‘W(i - ^Vr, * hr, * hhhr,’ 

Vr, - 


yj£„)/J= 


(A. 13) 


Vector Derivative Transformations 


Laplacian: 

vh . - 26fj^ + i,£^^)/j 2 + [tox^j - 2SXj^ + yV‘y5‘„ ‘ hU^ 

+ ( oy ^^ - 2 sy ^^ + yy „„)(«,£5 - «{£„) 1 / J ’ U . u ) 
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or 


■ (.fjj - 26t^^ + + tf^)/j2 (A.15) 

Gradient : 

?£ - Ky^fj - yjf^) i + (A.i6) 

Divergence: 

S ■ F . ly^(Fp^ - y^CFj)^ + Xj(rp_^ - x^(F,)jl/J (A.17) 

Curl: 

7 X F - k[y„(Fj)j - y^(Fj)^ - + x^(Fj)j]/J (A.18) 

Unit Tangent and Normal Vectors 
Normal to n-Line: 

n^^^ = (A. 19) 


Normal to ^"Line: 



n^^^ = V5/|VC| = (y i - X j)/»^ 
' H' n- 

(A.20) 

Tangent 

to n-Linc: 




t^*^^ E n<-^> 

xk - (x^i + y^4)/*^ 

(A.21) 

Tangent 

to 5-Line; 

t(^^ E n<^^ 

xk - -(x^i + 

(A. 22) 

Vector Components Tangent and Normal to C and n-Lines 



F E . 

n 

F = + x^F^)//y 

a.23) 


V (n) - (n) 

^t " ^ 

f “ ■*■ 

(A. 24) 


F E n^^^ • 

n 

F - (y_^Fj - x_F2)/^■ 

(A. 25) 


r (O * ,(0 
Ft - ? 


(A. 26) 
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Directional Derivatives 


3f/3n^’'^ s n^”^ • Vf - (yf - Bfp)/J»7 (A.27) 

- - n t 

3f/3t^^^ ^ t^’’^ • Vf - f^//7 (A.28) 

3f/3n^^^ = • Vf - (of. - Bf )/j/a (A.29) 

~ - ' 5 n 

3f/3t^^^ 2 t^^^ • Vf - - f /✓a (A. 30) 

- > - ri 
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APPENDIX B 


Concentration of Coordinate Lines Near a Body 

Consider the coordinate system generation Equations (2.3) applied 
to concentric circular boundaries of radius r^ and r^. With n ” 1 on 
the inner boundary, n ~ J on the outer boundary, and C varying mono- 
tonically from 1 to ^2 around these boundaries, the solution of 
Equations (2.3) can be given In the form 

X - r(n) cos [2 it (•^^-)] (B.la) 

y • r(n) sin I2ir ( 7 ^^-)] . (B.lb) 

Substitution of these expressions Into either of the equations of (2.3) 
with P(C,n) * 0 yields 

^ - ~ + Q(C.n) = 0 . (B.2) 

This can be made a perfect differential by taking the control function 
Q(C.n) of the form 

Q(C.n)H-|^ (B.3) 

where the minus sign has been introduced merely for convenience. 
Substituting Equation (B.3) into Equation (B.2) yields 


r 

r 


f" 

f 


which can be integrated twice to yield 

_(n) 


CzC 


Cj^f (n) 


(B.4) 


(B.5) 


The constants of integration may be evaluated from the boundary 
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conditions, r(l) * r^, r(J) - T2» so that 


r f(n) “ f(l) 

r(n) - {(r 2 /r^) . 

This equation nay then be solved for f(n) to yield 


(B.6) 




t(n) - f(l) - 
f(J) - f(l) r 


(B.7) 


ln[— ] 

>^1 

If the distance from the body to the Nth n-line is specified to be r^, 

the following equation must be satisfied: 

r 


f(N)__-_f(l) 


f(J) - f(l) r, 

Int— ] 


(B.8) 


It should be noted that Che form of f(n) is still arbitrary, subject to 


Equation (B.8). The form used in this study is given by 

f(n) - 

where K is a constant which must satisfy Che nonlinear equation 


(B.9) 


1 rl^ii 

N-l ^”^r 

NK ^ - 1 _ ^1 

1 1 r’^2, 

ln[— ] 


(B.IO) 


Once K has been determined from Equation (B.IO), the control function 
Q(5,n) is determined from Equations (B.3) and (B.9) to be 


■ - -(fr^ i”’' • 

*The form of f(n) was obtained from Reference [17). 


(B.ll) 


Although the above analyele was developed for circular boundaries, the 
effect will generally be the same for arbitrary boundary configurations. 

The procedure which gave the best results was to fix the distance 
of the n ■ 2 line at one one-hundredth of the boundary layer thickness 
from the body. Therefore, r^^ ■ r(2) is given by 


0.05 . 

r„ + r, . 


(B.12) 


It should be noted that r^ was set equal to the radius of a circle 
circumscribed about the body and that t 2 was set equal to the radius 
of the circular arc which represents a portion of the outer boundary. 
Now r^ and N have been specified so K is obtained by solving Equation 
(B.IO) using a false position iteration. 
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APPENDIX C 


r 


[ 

1 


( 

\ 




General Finite Difference Expressions 

This appendix suomarizes the finite difference expressions for a 
general function f(^,nit) which are necessary for implementation of the 
procedures outlined in Chapter IV. It should be noted that since 
the step size in the transformed plane le unity, it does not appear 
in the spatial difference expressions. For convenience, two short 
conventions are used in this appendix. All difference equations 
apply at the point denoted by the space subscripts (i,J) and the time 
superscript (n). The subscripts (1,J) and superscript (n) are 
understood where they are omitted. In this appendix, only time 
derivatives and C 'derivatives are shown. 


Time Derivative Approximations 
First-Order Backwards: 

f” “ (-f”“^ + f")/At 

Second-Order Backwards: 

f” * (f“‘^ - + 3f‘')/2At 

Spatial Derivative Approximations (Second-Order Only) 
First Derivative, 

First Derivative, 


Central Difference: 


“ ^i+l,j ■ ^i-l,j 
Forward Difference: 






(C.2) 


(C.3) 


(C.A) 
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First Dsrlvstiv*, Backward Dlffaranca: 

" ^^1-1 ,J 

Second Derlvaciva, Central Differences 

‘ ^1+1, j ■ 


(C.5) 


(C.6) 


^Kr\ “ ^^i+l,J+l " ^1+1, J-1 ■ ^i-l,J+l ^1-1, J-1' 


Cross Derivative, Central Difference: 

,)/4 (C.7) 

Cross Derivative, C-Forward Difference and n-Central Difference: 

■ *"*i+2,J+l * ^l+2,j-l ***i+l,J+l ■ *i+l,3-l) 

Cross Derivative, C-Backward Difference and n*Central Difference: 
^*l-2.J+l ■ *1-2,J-1 ■ ***1-1,J+1 ■ 

Cross Derivative, C-Central Derivative and n-Forward Difference: 

f s (f - f ^ L(f • f 

' i-l,j+2 ^i+l,j+2 ^ ^*^i+l,j+l ^i-l,j+l) 


”^^^i+l,J “ ^1-1, 


(C.IO) 
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